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The first results on lattice at finite temperature T and chemical potential fj, with 2+1 flavors 
of Wilson quarks are presented. We have studied the dependence of chiral phase transition and 
deconfinement phase transition on quark mass. Finite volume size analysis and Binder cumulants 
are used to determine the properties of phase transition. Phase diagram of QCD with 2+1 flavors 
of Wilson quarks are presented. 
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I. INTRODUCTION 

With the heavy-ion collision experiments at RHIC now 
running, it is much more important to understand the 
phase diagram of QCD at finite temperature and fi- 
nite density. QCD predicts the existence of quark-gluon 
plasma(QGP) at high temperature and the existence of 
phase transition between the hadronic matter and the 
QGP. But there is no direct evidence yet. It is a very 
important work to locate the phase transition region and 
to investigate its properties. 

The lattice gauge theory (LGT) is the unique nonper- 
turbative method for the quantum gauge field. Lots of 
MC (Monte carlo) simulations for QCD on lattice were 
performed to study the quark mass spectrum, gluon ball 
mass spectrum, phase diagram and so on. 

In the past, the most MC simulations for QCD at fi- 
nite temperature and finite density were performed with 
staggered quarks owing to its preserving the continuous 
U(l) chiral symmetry . But staggered fermions formal- 
ism on the lattice does not completely solved the species 
doubling problem. For to reduce the fermionic degrees 
of freedom, MC simulation take the fourth root of the 
fermion determinant to thin the species doubling prob- 
lem. The operation can arise the locality problem in 
simulations with non 4 taste fermions The Wilson 
fermion formalism preserves the flavor symmetry and has 
no the species doubling problem. It is also only known 
as a formalism that possesses a local action for any num- 
ber of flavors and is well formalism for MC simulation on 
lattice. 

In our previous workQ, the first phase diagram of 
QCD with 4 flavors(./Vf = 4) of Wilson quarks at small 
chemical potential and high temperature was presented 
by using imaginary chemical potential method. The crit- 



ical line of first order phase transition was located at 
(/i, T) plane. Nf — 4 presume that there are four de- 
generate quarks. But the physically relevant case should 
be two degenerate massless u,d quarks and one light S 
quark, which is so-called 2+1 flavors quarks. It is im- 
portant to investigate the phase diagram of QCD with 
2+1 flavor quarks for us to understand the actual state 
of high density star and the cosmic origin. 

In this paper, we will present our results of MC simula- 
tions with 2+1 flavors of Wilson quarks at finite chemical 
potential and finite temperature by using the pure imag- 
inary chemical potential method. 

In section II, we will briefly instruct the latterly 
progress in this domain. In section III, we will presented 
the formalism used in this paper. In section IV, our de- 
tailed results of MC simulation were presented. At last, 
we will present the conclusion. 



II. GENERAL CHARACTERISTIC OF THE 
PHASE TRANSITION FOR QCD 
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FIG. 1: The schematic phase diagram of QCD in (/x, T) plane 
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Before starting our discussion, let us briefly review 
present understanding of the phase diagram structure. 
The phase diagram in the (//, T) Plane is shown in FigQ] 
The (p,T) Plane is divided into three main parts, I, II, 
III. At I domain, where the temperature and chemical 
potential are below critical line, the QCD matter was 
formed as hadrons and the quarks are confined. At the 
II domain, because of extra high temperature and strong 
interacting, quarks are released and the strongly interact- 
ing QCD matter formed the quark-gluon plasma(QGP) 
and the quarks are deconfined. In the III domain, where 
the temperature is low and the density of quarks is very 
high, it is believe that two quarks on the fermion surface 
can be bound into an Cooper pair, which is so-called color 
superconductivity. That case is similar to the supercon- 
ductivity in the electromagnetism. The solid line which 
end at the point A is the the first order critical line. The 
dotted line which start from the point A is the the second 
order critical line. The point A is the trecritical point. 
It is very important to accurately locate the trecritical 
point for searching and investigating the QGP. 
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FIG. 2: The schematic phase diagram of QCD in (m u d,m 3 ) 
plane for /i = 

According to the MC simulation result of the past on 
lattice, the nature of phase transition of QCD depends on 
the quantity of quark mass. At zero chemical potential 
(i = 0, the schematic phase diagram is shown in FigJ2]in 
the {m u d,m a ) plane. For the pure gauge case, for which 
all the quark mass degenerate into infinitely large, a first 
order transition is expected in the quenched lattice gauge 
theory. For QCD with two degenerated light u,d quarks 
and a infinitely large s quark, which is so-called two fla- 
vors QCD, it is suggested that the chiral transition is 
of second order in the limit of vanishing quark massQ. 
For three or more degenerate quarks, the transition is 
characterized as first order while the quantity of quark 
mass is less than some critical value m c q Q ■ For the non- 
degenerate quarks, the strange quark is heavier than the 
u,d quark, one expects a critical line which separates the 



regime of first order phase transitions from the crossover 
regime pi- In the other area, the transition is continuous. 
The structure of phase diagram is sensitive to the the 
non-degenerated quark mass. 

III. THE FORMALISM IN LATTICE GAUGE 
THEORY 

In the thermal field theory, the canonical partition 
function Z(V, T) can be expressed as: 

Z = J [dU][d^][dij:]e~ ( ~ SG ^ +s ^ u ^^\ (1) 

where Sq is the gauge field action: 

1 f 1/T f 
S G=zJ dT J dZ * F ^ F w 

Sf is the fermion field action: 
f 1/T f " 

JO J v j 

where m/ is the quark mass. 

On lattice, the canonical partition function Z is: 

f Nf 

Z= [dU] Y[ detM(U,K } ,an)e- pSG[u \ (2) 
f= i 

where, k/ = 1/(8 + 2am f) is the quark mass parameter, 
a is the lattice space, /i is chemical potential. The gauge 
field action Sq is: 

s G (u) = -pY1 p w ( 3 ) 

where P M)I , is a plaquette on lattice. For QCD with Wil- 
son fermions, the fermion action Sf is : 

S F ($,il>,U) = $Mil>, (4) 

where M is the fermionic determinant: 

3 

Mij = 5 i)i +K / {^[(l- 7 i)C/» + (l + 7i)^(ft-«)] 

i=l 

+ [(1 - l4 )e a »U 4 (n) + (1 + 7i )e-^Ul(h - 4)]}. 

(5) 

For /i = 0, detM(U, Kf, a/x) is positive real, MC methods 
can be used to simulate the thermodynamics of QCD. 
But for /i ^= 0, detM{U, Kf, afj,) is complex, the standard 
MC simulation methods are invalid. Some improved MC 
methods can be introduced in simulation of QCD . In this 
paper, the pure imaginary chemical potential method is 
chosen, which is proper for low baryon density and high 
temperature. 
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QCD is a theory of asymptotic freedom. The strength 
of interaction between two static quarks could reflects 
on the state of quarks. Strong interaction indicate that 
quarks are isolated and deconfined. Polyakov loop 



P{x) = Tr 



N t -l 

n ut&t) 



(6) 



can tell us the strength of interaction between two static 
quarks, which can be used to determine the confinement - 
deconfincment phase transition. The chiral condensate 



1 



[dU}[dip][di>]^e 



-S„ — S r 



1 



ZVN t 



[dU}Ti (M-^U}) (DetM[U]) N ' e 



(7) 



is the order parameter of chiral-symmetry breaking. At 
the chiral limit, non-zero of (ipip) implied that the chiral- 
symmetry of QCD is spontaneously broken. According 
to the present understanding, the chiral transition point 
is same as the deconfinement transition one. 

In this paper, several methods are used to locate the 
critical point of transition. A method is to search the lo- 
cation of rapid variety of order parameters. The second 
one is to search the peaks of the order parameter fluctua- 
tions at critical points. It is convenient to determine the 
universality class of transition by using the Binder cu- 
mulants of physical quantity at the transition point. The 
dimensionless fourth order cumulants B4 of an arbitrary 
observable physical quantity x is constructed as 



-B 4 (m c ,^ c ,/3 C ) 



< {Sx) 4 > 

< (5x) 2 > 2 



(8) 



where, m c ,fj, c ,/3 c is the value at critical point. Binder 
cumulants calculated on different size lattices will inter- 
sect at some value of the quark mass. At the infinite 
volume limit, these intersection points will converge to a 
universal value which is characteristic for the universality 
class. The universal value is quite different for different 
symmetric system. In Nf — 3, the chiral symmetry of 
QCD at critical point is like 3D Ising model, the inter- 
section points converge to B 4 = 1.604 @. In this paper, 
-B4 is used to determine the characteristic of chiral phase 
transition. 



IV. SIMULATION RESULTS FOR QCD WITH 
2+1 FLAVORS FERMIONS 

In this paper, MC simulation are performed with the 
standard R algorithm^. We have modified the MILC 
collaboration's public codeQ to fit the imaginary chem- 
ical potential and 2+1 flavor case. MC simulations are 
performed at lattice size V x N t = 8 3 x 4 and finite size 
scaling analysis was checked for some different lattice size 



in certain parameters. The parameter ki, which corre- 
sponds to parameter of u, d quark mass, was scanned 
from 0.16 to 0.20. The parameter K2, which corresponds 
the mass parameter of s quark and is heavier than u, d 
quark, was fixed at k,i — 0.08 which is much less than 
K\. The St is 0.02. Each configuration runs 20 molec- 
ular steps with microcanonical step size St — 0.02. A 
configuration is token for measurement every other 20 
configurations. 
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FIG. 3: Polyakov loop norm as a function of for ajii 
0.0,0.15,0.20,0.24 in k x = 0.16. 




FIG. 4: Susceptibility of Polyakov loop norm as a function of 
P for afii = 0.0, 0.15, 0.20, 0.24 in ki = 0.16. 

At First, we fix u,d quark mass parameter K\ at 
0.16, which is far from chiral limit. MC simulations 
are done by scanning the parameter /3 for different a/ij. 
2000 ~ 4000 configurations are measured for every sim- 
ulation. Polyakov loop {\p(x)\) and its susceptibility 
X|p(s)| as function of (3 are shown in Figl3] and Figj4] 
As one can see in FiglHl with (3 increased, (|p(x)|) raise 
from approximate zero to obvious non-zero. Especially at 
j3 ~ 5.25 ~ 5.35, the Polyakov loop show rapidly growth. 
Exactly at the same location, X|p(a;)| show sharp peaks for 
different a/x/ in Fig(5] The points at the peaks of X|p(x)| 
should be the critical points of deconfinement phase tran- 
sition. At the same condition, the chiral condensate was 
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also studied synchronously. The chiral condensate (i/jip) 
and the susceptibility x^ are illustrated in Fig[5] and 
FigJS] As P is increased, show decrease in FigJSJ At 
the peaks of X\p(x)\ hi Fig|4] Xm show rapid drop and the 
susceptibility show sharp peaks too. The results implied 
that the deconhnement phase transition maybe coincides 
with the chiral phase transition. On the other hand, it 
is obvious in Figj5] that the chiral condensate ("ipip) is 
nonzero while QCD is in the deconhnement phase. It 
is due to non-preserving the chiral symmetry for Wilson 
fermion action . That means the chiral condensate is not 
any more a good order parameter for the chiral transition 
of QCD with Wilson fermions. However, we can also find 
the chiral critical point by locating its rapid change and 
the sharp peak of its fluctuation. 



term: 




FIG. 5: Chiral condensate as a function of /3 for a/ij 
0.0, 0.15, 0.20, 0.24 in /ti = 0.16 and k 2 = 0.08. 




FIG. 6: Susceptibility of chiral condensate as a function of /3 
for afn = 0.0, 0.15, 0.20, 0.24 in Ki = 0.16 and k 2 = 0.08. 

In Fig[J]and Fig[51 the critical points for different a/.ij 
can be located approximately. All critical points can be 
collected into a critical line {Pc,a^i) m ( a M/j/3) plane. 
For small imaginary chemical potential, this critical line 
could be expand into Taylor series. From our simulation 
results, we obtain the critical line up to the quadratic 



(3 C = 5.264(1) + 0.75(4) (a^) 2 + O (a 4 fij) 



(9) 



The fitted line is shown in FigJTl By replacing /i/ with 
—i/j., the critical line © is transformed from imaginary 
chemical potential case to real chemical potential case: 



P c = 5.264(1) - 0.75(4) (a/i) 2 + O (aV) 



(10) 




FIG. 7: The fitted deconhnement phase transition line by 
collecting the critical points 

According to the renormalization group relation be- 
tween the lattice spacing a and P, the two loop pertur- 
bative expression gives 



a Al = cxp 



4tt 2 



33 - 2Nf 
459 - 57A/ 



P 



(11) 



In 



8n 2 



-P 



(33-2A/) 2 \33-2N f 



where Al is the lattice QCD scale. The temperature T 
is related to a and Nt by T — \/(aNt). The critical line 
on the (fj,, T) plane is shown in FigO 

In order to determine the properties of the transition, 
the dependencies of the transition on the spatial volume 
of the lattice are studied. While the lattice volume is 6 3 x 
4, 8 3 x 4, 12 3 x 4, respectively and m is fixed at 0.16, the 
results are shown in FigJ^l It is clear that the fluctuation 
°f X\p(x)\ have no obvious dependence on the the lattice 
size and the location of critical points have not shown 
distinguishable shift. It suggests that the transition is 
just crossover while the u,d quark mass is far away from 
the chiral limit. 

At the following, the dependences of phase transition 
on the quantity of u,d quark mass near the chiral limit 
are studied. Here, the chemical potential a\ii is fixed 
at 0.15, we scanned u,d quark mass parameter k\ from 
0.16 to 0.180. As K\ is increased, which correspond to 
decrease the quantity of u, d quark mass, the results of 
MC simulation at the critical points show likewise change 
as above paragraphs. The detailed results of (\p(x)\) and 
X\p(x)\ versus P were shown in FigfiTJl and FigfTTI 
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FIG. 8: Phase diagram on the (fiB,T) plane. The area be- 
tween the two dashed lines is the error band for 4 flavor Stag- 
gered quarks at m, = 0.05 fromQ], the dotted lines is for 
2 flavor Staggered quarks at m q = 0.025 from [J, and the 
solid lines is for 2+1 flavor Wilson quarks at ki = 0.16 and 

K 2 = 0.08. 
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FIG. 9: Susceptibility of Polyakov loop norm as a function 
of f3 for different lattice volume V = 8 3 ,I2 3 ,I6 3 in Ki = 
0.16, k 2 = 0.08. 




FIG. 10: Polyakov loop norm as a function of j3 for different 
fti in ajii = 0.15 




FIG. 11: Susceptibility of Polyakov loop norm as a function 
of j3 for different ki in a\ii = 0.15. 

It is obvious in Fig[l0l for large value of K\ the as- 
cending rate of the quantity of (\p(x)\) near the critical 
points are more steep and the peaks of X\p(x)\ m Fig llll 
at the critical points are more sharp and higher. That 
implies that for tiny u,d quark mass the deconfincmcnt 
phase transition could be of first order. 

The finite volume size effects are checked in K\ = 0.170 
and afii = 0.150. 4000 ~ 10000 configurations are mea- 
sured for every simulation. The Xb(x)| as function of 
(3 are shown in FigJT2] The peak value of X\p(x)\ ^ s 
shown higher for greater lattice volume and the location 
of the critical point is shown obvious shift with differ- 
ent lattice volume size. For a second order transition, 
Xmax oc V a , a < 1. Our results are fitted to Xmax oa V a , 
and we get a w 0.5. The result suggests that the transi- 
tion is of second order in the situation. 




FIG. 12: Susceptibility of Polyakov loop norm as a function of 
(3 for lattice volume V = 8 3 , 12 3 , 16 3 in m = 0.170, k 2 = 0.08 
and a/j,i — 0.150. 

As one can see in FigOJJ for Ki > 0.175 the variety of 
(\p(x)\) at the critical points is very rapid. In k\ = 0.175 
and (3 = 5.10, the history of MC simulation is shown 
in FigQ2] and the histogram of the Polyakov loop norm 
is shown in Fig03] Two plateaus at \p{x)\ ~ 0.25 and 
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\p(x)\ « 0.55 alternately appeared in the history of MC 
simulation and two peaks respectively at \p(x)\ « 0.25 
and \p(x)\ « 0.55 are distinguishable in FigQT] Espe- 
cially for tiny u, d quark mass , the structure of two peaks 
in the histogram are shown more outstanding. The re- 
sults suggest, while the quantity of u, d quark mass is 
decreased to be close to the chiral limit and the s quark 
mass is large enough, the deconfinement phase transition 
should be of first order. 
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FIG. 13: History of MC simulation for Polyakov loop norm 
in afj.i — 0.15, k1 — 0.175 and K2 = 0.08 
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FIG. 14: Histogram of Polyakov loop in a{ii = 0.15, k1 — 
0.175 and K2 = 0.08, the solid curves are fitted Gaussian 
curve. 

In order to investigate the universality class of the 
phase transition, Binder cumulants are calculated in dif- 
ferent quark mass and chemical potential. For K2 = 0.08, 
MC simulation are performed by scanning am- 2000 ~ 
8000 configurations are measured for every simulation. 
The results of B4 are shown in FigfTSJ With increas- 
ing am, B4 drops. For am — * 0, BA « 2.0. A sec- 
ond order point, which is chosen as the critical value 
BA = 1.604 of 3D Ising model, separate the first order 
from the crossover. The phase transition for BA < 1.604 
should be of first order and for BA > 1.604 could be 



crossover. According to above results, the phase tran- 
sition near BA w 1.604 should be of second order and 
for BA ^> 1.604 should be crossover, their exact border 
can not be located by traditional methods. From FigfT5l 
we got the second order critical cure in Fig |16l in the 
(a/i/jKi) plane. For small a/xj, we can Taylor-expand 
the critical line in Fig[T|5] as: 

Ki{am) = K i i a m = 0) + b(am) 2 + 0{a 4 m 4 )- 

From FigOS we get k?(0) = 0.1864(4), b = -0.280(3). 
By replacing m by — ip,, the critical line is continued from 
imaginary chemical potential to real chemical potential: 

K c x {am = 0.1864(4) + 0.280(3)(a^) 2 + O(aV)- (12) 

According to (fT5|) . for small n-y (large u,d quark mass) 
the critical point A in FigQ] is more close to the T axis. 
Especially while K\ = /ti(0), the critical point A will 
collide with T axis and disappear. The result coincides 
with 0. 




FIG. 15: The Binder cumulant B4 of Polyakov loop norm 
for different (a(j,i) 2 . The dashed line corresponds the critical 
line B4 = 1.604. The dash_dotted line are the fitted lines for 
different ki 
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FIG. 16: The second order critical points in the (a/i/,Ki) 
plane. 
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V. CONCLUSIONS 

In this paper, we have used MC simulation on lattice 
with 2+1 flavors Wilson fermions to determine the struc- 
ture of the phase diagram of QCD at finite temperature 
and small density. As we have discussed in the part of 
instruction, the operation that the fermionic determinant 
is token the fourth root to thin the species doubling prob- 
lem in MC simulation with staggered fermions can arise 
the locality problem in simulations. Whereas the MC 
simulations with Wilson fermions on lattice have no this 
puzzled problem. 

According to our results of MC simulation, the 
schematic expected phase diagram was shown in Fig |17l 
Below the chiral limit surface k u ^ = n c hirai and above the 
surface K u d = 0, the space was separated into two parts 
by the phase transition interface. QCD is in hadronic 
matter phase inside the phase transition interface and in 
QGP phase outside the phase transition interface. The 
phase transition interface is separate into the first order 
and the crossover by a second order critical curve. The 
critical curve is just the critical point A in Fig[l] which is 
function of quark mass. 



FIG. 17: The schematic phase diagram of QCD with 2+1 
flavors of Wilson quarks in the (fi,T,K u d) space. A critical 
line separate the phase transition interface into the first order 
and the crossover. The phase transition is of second order 
near the critical line. 
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